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Abstract The purpose of this paper is to construct the early exercise boundary for 
a class of nonlinear Black-Scholes equations with a volatility function depending 
on the option price. We review and revisit a method how to transform the problem 
into a solution of a time depending nonlinear parabolic equation defined on a fixed 
domain. An example of numerical computation of the early exercise boundary for a 
nonlinear Black-Scholes equation is also presented. 



1 Black-Scholes equations with a nonlinear volatility function 

The main purpose of this paper is to review and revisit the fixed domain transfor- 
mation method adopted for solving a class of nonlinear Black-Scholes equations 
having the form: 



dt 



+ { r - q )S— + -o z {S l d l s V 1 S,T-t)S l j- I - r V = 0, S>0,te(0,T). (1) 



A solution V = V(S,t) can be identified with a price V of the option contract in the 
future maturity time T > (e.g. call or put) where S > is the underlying asset 
value at the present time t G [0, T). Here, r > is the riskless interest rate, q > is 
the dividend yield rate of the underlying asset. For American style of a call option, 
the free boundary problem consists in construction of the early exercise position 
Sf = Sf(t) and the solution V = V (5, t) to equation (Q~|i defined on the time dependent 
domain < S < Sf(t), < t < T (cf. Kwok 02)). V is subjected to the boundary 
conditions yielding C 1 smooth pasting of V (S, t) and V(S,T) at S = Sf(t): 
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7(0,0 = 0, V(S f (t),t)=S f (t)-E, d s V(S f {t),t) = l, (2) 
and the terminal pay-off condition at expiry t = T, 

V(S,T) = (S-E) + , (3) 

where E > is the exercise price. 

We briefly mention a motivation for studying the nonlinear Black-Scholes equa- 
tion having the form of (fTJ. Such equations with a volatility function a(S 2 dgV, S,T — 
t) arise from modeling the option prices by taking into account nontrivial transac- 
tion costs (cf. Leland lfl4l . Hoggard et al. iflOl . Avellaneda and Paras Q), market 
feedbacks and effects due to large traders choosing given stock-trading strategies 
(Frey Q, Frey and Patie (8), Frey and Stremme [9], During et al. I0, Schonbucher 
and Wilmott |[T5l ), the risk adjusted pricing methodology model due to Kratka lfl2l 
and its modification developed by Jandacka and Sevcovic 11 111171 ). As an example 
for application of the numerical method, we consider on a nonlinear model taking 
into account imperfect replication and investor's preferences which has been pro- 
posed by Barles and Soner in J4). If investor's preferences are characterized by an 
exponential utility function they derived a nonlinear Black-Scholes equation with 
the volatility function a given by 

G 2 (S 2 d 2 V,S, t) = a 2 (1 + W{a 2 e n S 2 d 2 V)) . (4) 

Here a 2 > is a constant historical volatility of the asset price returns, f is the 
unique solution to the ODE: f" (x) = (W(x) + 1 )/ (2 - *) , f (0) = and a > 

is a constant depending transaction costs and investor's risk aversion parameter 
(see H for details). The function f satisfies: f(x) = 0(x%) forx and f (x) = 
0(x) for x — > oo. For practical purposes, the solution f (x) can be constructed from 
an implicit equation obtained in . 

We revisit an iterative numerical algorithm for solving the free boundary prob- 
lem ([TJ-ffSJi developed by Sevcovic in |17|. The key idea of this method consists 
in transformation of the free boundary problem into a semilinear parabolic equa- 
tion defined on a fixed spatial domain coupled with a nonlocal algebraic constraint 
equation for the free boundary position. This method has been analyzed and utilized 
in a series of papers lH~ll2l lT6l - |T9i by Ehrhardt and Ankudinova and the author. The 
disadvantage of the original method consists in the necessity of solving an algebraic 
constraint equation. In this approach, highly accurate evaluation of the derivative of 
a solution at one point entering the algebraic constraint is needed (cf. ifPTl ). In this 
note, we present a new efficient way how to overcome this difficulty by considering 
an equivalent integrated form of the algebraic constraint. We also present results of 
numerical calculation of the free boundary position for the Barles and Soner non- 
linear extension of the Black-Scholes model. 
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2 Fixed domain transformation of the free boundary problem 



We recall the method how to transform the free boundary problem (Q])-© into a 
form of a nonlinear parabolic equation defined on a fixed domain and satisfying a 
nonlocal algebraic constraint equation developed by the author in ifTTl . It is based 
on the following change of independent variables and the transformed function IT = 
n(x, z) defined as follows: 

r = T-t, x = ln(p(z)/S), n(x,z) =V(S,t)-Sd s V(S,t), (5) 

where p(z) = S f (T - z). Clearly, T 6 (0,T) and x E (0,°°) iff S E (0,5/0)). The 
boundary value x — corresponds to the free boundary position S = Sf(t) whereas 
x = +°° corresponds to the default value S = of the underlying asset. Under the 
structural assumption 

0<q<r 

made on the interest and dividend yield rates and following derivation of the equa- 
tion for n, it turns out that the function II and the free boundary position p satisfy 
the following system of parabolic equation © with algebraic constraint @: 

dn ( , s a 2 \dn id/ 2 dn\ 



n(o, t) = -e, n(+oo, T ) = 0, jc> o , z e (0,r) , 

-E foix < ln(r/q), 
0, otherwise, 



n(x,o) = 



g + ^ g ^ M M with p(0)^, ,7, 

q Zq ox q 

where a 2 = a 2 (d x n(x 7 z),p(z)e- x ,z) , = + r-« (cf. Q3). Notice that 
equation (|7]) is not quite appropriate for construction of a robust numerical approxi- 
mation scheme since any small inaccuracy in approximation of the value <9 v n(0, t) 
is immediately transferred in to the entire computational domain x € (0,°°) through 
the free boundary function p (t) entering ©. Instead of (0, we present a new equiv- 
alent integrated equation for the free boundary position p(z). Indeed, integrating the 
governing equation (O forx S (0,°°) taking into account the boundary conditions 
n(Q, z) = —E,II(°°, z) = (and consequently <9 v n(°°, z) = 0), we obtain the fol- 
lowing spatially integrated form of the algebraic constraint: 

1 (Elnp(z)+ ( II(x,z)dx) +qp(z)-qE 



dx \ Jq 



/o 



- X - a 1 (d x n(x, x) lP {x)e-\x)^- (x, z) + rn{x, z) 
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3 Numerical scheme based on operator splitting technique 



The idea of the iterative numerical algorithm is based on the original numerical 
discretization scheme proposed by the author in [17]. We modify this method by 
considering the alternative integrated form ^ of the constraint between TI and p . 
The spatial domain x E (0,°°) is restricted to a finite interval of values x E (0,L) 
where L > is sufficiently large. For practical purposes one can take L w 3 (see 
|[T7lO . Let us denote by k > the time step, k = T jm and by h > the spatial 
step, h = L/n where m,n G N stand for the number of time and space discretization 
steps, respectively. We denote by U- an approximation of n(x,,T/), p 7 ~ p(T;), 
b J K,b(Xj) where x, = ih, X\ = jk. We furthermore denote by TV the vector W = 
{TI- ,i = 1, . . . ,«}. We approximate the value of the volatility a at the node (xi, Tj) 
by the finite difference approximation as follows: 

We set n?(x) = TI(xi,0). Next, following the idea of the operator splitting method 
discussed in [17), we decompose the above problem into two parts - a convection 
part and a diffusive part by introducing an auxiliary intermediate step 17 J_ 2. Our 
discretization of equations ^ and (O reads as follows: 
(Integrated form of the algebraic part) 

E In p i = E In p j ~ 1 + 7 {Tl^ 1 ) - I (W )+k(qE- qp j - h (p j , n j ) ) , (9) 

where Io(TI) stands for numerical trapezoid quadrature of the integral /J° TI(%)d% 
whereas I\ (p 7 ,i7) is a trapezoid quadrature of the second integral in (8j, i.e. 

poo / 2 dTI \ 

Ii(p J ,n)^ J (--a 2 (d x n(x),p J e- x : Z j ) — (x)+rn(x)\d x . 



(Convective part) 



(Diffusive part) 



+b J —n J -i=o, (io) 

k ox 



The convective part can be approximated by an explicit solution to the transport 

y 1 

equation d T TI + b(x)d x TI = 0. Thus the spatial approximation iZ 1 can be con- 
structed from the formula 

n J-k = in j - 1 m if^=x i -lnpJ+ln P J- 1 -(r- q )k>0 1 
' 1 —E, otherwise, 
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where a piecewise linear interpolation between discrete values Tlj _1 ,z = 0, 1, ...,«, 
is being used to compute the value IT- 7-1 (xj — lnp 7 + lnp 7-1 — (r — q)k). 

The diffusive part can be solved numerically by means of finite differences. Using 
central finite difference approximation of the derivative d x W we obtain 



k 2 2/i 

Now, equations (O, dT2l > and ( foi l can be rewritten in the operator form: 

p j = &(n j ,p j ), n j -i = &(n J ,p j ), ^(n j ,p J yn j = n j -i, 

where &(JV >p J ) is the right-hand side of the integrated algebraic equation (O. 
The operator ^{Tl J ,pi) is the transport equation solver given by the right-hand 
side of (fT2l and srf = srf{H ] \p J ) is a tridiagonal matrix with coefficients given 
corresponding to H3i . At each time level tj,j = l,...,m, the above system can be 
solved approximately by means of successive iterations procedure. Given a discrete 
solution IP -1 , we start up iterations by defining IP' = fl ;_1 ,p ; = p-' -1 . Then 
the (p + l)-th approximation of fl 7 and p ' is obtained as a solution to the system: 

p7>+i = &(W'P,pi>P), W-^ p+l = ^(W' p ,p j - p+l ), 

£/(W'P,pj'P +1 )W'P +1 = n J -*' p+1 . (14) 

We repeat the procedure for p = 0,\...,p max , until the prescribed tolerance is 
achieved. 

At the end of this section, we present a numerical example of approximation of 
the early exercise boundary for the Barles and Soner model by means of a solu- 
tion to the transformed system of equations. In this model the volatility is given by 
expression (0). A discrete solution pair (iT,p) has been computed by our iterative 
algorithm for the model parameters: E = 10, T = 1 (one year), r = 0.1 (10% p. a) , 
q = 0.05 (5% p. a.) and a = 0.2. As for the numerical parameters, we chose n = 750 
spatial points and m = 225000 time discretization steps. The step k = T jm repre- 
sents 140 seconds in the real time scale. In order to achieve the precision 10~ 7 we 
used Pmax = 6 micro-iterates in ( fT4b . A graphical plot of the early exercise bound- 
ary p (t) = Sf(T — t) is shown in Fig.Q] Taking a positive value of the risk aversion 
coefficient a = 0.15 resulted in a substantial increase of the free boundary position 
p(r) in comparison to the linear Black-Scholes equation with constant volatility 
a = a. Notice that the Barles and Soner model for a = coincides with the linear 
Black-Scholes model with constant volatility. 
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Fig. 1 A comparison of 
p(T) = Sf(T - t) (solid line) 
for the Barles and Soner 
model with a = 0. 15 and for 
the Black-Scholes equation, 
i.e. a = 0. 
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